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We propose a simple extension of the well known ST2 model for water [F.H. Stillinger and A. 
Rahman, J. Chem. Phys. 60, 1545 (1974)] that allows for a continuous modification of the hydrogen 
bond angular flexibility. We show that the bond flexibility affects the relative thermodynamic 
stability of the liquid and of the hexagonal (or cubic) ice. On increasing flexibility, the liquid-liquid 
critical point, which in the original ST2 model is located in the no-man’s land (i. e. the region where 
ice is the thermodynamically stable phase) progressively moves to a temperature where the liquid 
is more stable than ice. Our study definitively proves that the liquid-liquid transition in ST2 is a 
genuine phenomenon, of high relevance in all tetrahedral network-forming liquids, including water. 


The possibility that a one-component system assumes 
(beside the gas phase) more than one disordered con¬ 
densed phase is currently highly debated in liquid state 
physics. Since the original proposition [Ij of a liquid- 
liquid (LL) transition in supercooled water (based on a 
molecular dynamics study of the ST2 Q potential), a 
large literature body has investigated this topic, suggest¬ 
ing that the microscopic origin of a LL transition must 
be attributed to the competition between two local struc¬ 
tures, differing in energy, entropy and density [3-14|. 
Still, when and how the interaction potential between 
molecules will favor a LL transition which can be ac¬ 
cessed in the absence of crystallisation is rather unclear. 
Only recently an effort has been made to provide a pic¬ 
ture which simultaneously accounts for the free energy of 
both the liquid and ordered phases as well as of 

the kinetic barrier separating them. The driving force be¬ 
hind this renewed interest in the physics of LL transitions 
has been provided by two very controversial studies from 
the same group 15|,ll8|- These studies state that in previ¬ 


ous numerical investigations — including the ST2-based 
results which originated the LL transition concept — the 
low density liquid phase appearing below the LL criti¬ 
cal point was in reality an ice phase, i.e. crystallization 
was mistaken for a LL transition. Several following in¬ 
vestigations by different groups have disagreed with this 
interpretation, providing further support in favour of the 
presence of two well-d efined distinct liquid phases in the 
ST2 model 


The most recent contribution [Ifilj has confirmed that 
the free energy basins of the two liquids are well sepa¬ 
rated from the crystal one and hence, in principle, both 
liquid phases can be explored in metastable equilibrium 
conditions. Of course, this requires that the metastable 
liquid phases survive for a time longer than the equili¬ 
bration time. Such times can not be calculated by ther¬ 


modynamic information only. It is thus worth exploring 
the possibility of a definitive proof of the existence of a 
LL critical point in ST2 that does not require kinetic in¬ 
formation. We present such a proof here by continuously 
tuning one of the ST2 model parameters. We show that 
it is possible to modulate the relative stability of the liq¬ 
uid and of the hexagonal (or cubic) ice ly/c such that the 
melting temperature of T/c drops below the LL critical 
temperature. Under these conditions, the low-density liq¬ 
uid is thermodynamically more stable than I^/c, demon¬ 
strating that these two phases are distinct. The results 
reported in this Letter not only conclude once and for all 
the debate on the existence of a genuine LL transition 
in the ST2 model but also provide important clues on 
the mechanisms controlling crystal stability in tetrahe¬ 
dral lattices. 


Our study builds on recent investigations of patchy 
colloidal particles, interacting via four attractive patches 


tetrahedrally located on the particle surface III 
Searching for the particle properties favoring the self- 
assembly of the technologically relevant diamond lat¬ 
tice [ 23 , it has been discovered that very flexible bonds 


destabilize open crystal phases so much that the liq¬ 
uid retains its thermodynamic stability even at very 
low temperatures [ 2 ^. Crystallization is instead favored 
by highly directional bonds. In addition, at low tem¬ 
perature T these systems can exhibit a LL transition 
between two interpenetrating tetrahedrally coordinated 
networks 17|. On increasing the bond flexibility the 
LL transition becomes thermodynamically stable. These 
results have been confirmed in another colloidal model 
mimicking DNA constructs with valence four 26|. Both 
colloidal models are characterized by a very large inter¬ 
particle softness, allowing for full network interpenetra¬ 
tion [23 ■ As a result, the density of the coexisting high 
density liquid approximately doubles the density of the 
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low density liquid phase, a factor significantly larger than 
what is expected for water. This extreme softness casts 
some doubt on the applicability of these results to molec¬ 
ular systems and water in particular. We alleviate these 
doubts here, supporting once more the hypothesis that 
the liquid-liquid transition is a genuine feature of tetra¬ 
hedral network forming liquids. 

The original ST2 potential envisions a water molecule 
as a rigid body: the oxygen atom (O) is located at the 
centre, while the two protons (H) are located at a dis¬ 
tance of 1 A from O, forming a fixed HOH tetrahedral an¬ 
gle. Two sites (mimicking the lone-pairs, LP) are located 
at distance 0.8 A from O, such that the two 0-H and 
the two 0-LP unit vectors form the vertices of a perfect 
tetrahedron. The H and LP sites carry an electric charge. 
Long-range electrostatic interactions are included via the 
reaction-field. Complete details of the simulation proce¬ 
dure are as described in Ref. [![. For this model, the 
phase diagram has recently been evaluated, demonstrat¬ 
ing the stability of ice \hjc at low temperature and pres¬ 
sure [ 2 ^ , consistent with the recent observation of ice I in 
simulation of ST2 [ 2 ^. We modulate the flexibility of the 
hydrogen bonds by allowing the unit vectors pointing to¬ 
ward the H and LP sites to fluctuate (with no additional 
energy cost) with respect to the original direction, with a 
maximum angle 0max (see Fig. [T]). By changing ^max it is 
possible to continuously tune the bond flexibility. When 
^*max = 0° the modified model coincides with the original 
ST2 model. Apart from 0„iax = 0° (cos^max = 1-0), we 
explore in detail the cases ^max = 8.11° (cos^max = 0.99), 
^*max = 11.5° (cos^max = 0.98) and the case 0i„ax = 14° 
(cos^max = 0.97). We note that in principle, an energy 
cost to bending could be included in the model. How¬ 
ever, the main effect of this would be to make the effec¬ 
tive bond flexibility temperature-dependent, and we have 
thus omitted this here. 

To provide evidence that on increasing 0max, the 
tetrahedral network becomes more and more flexible we 
evaluate the 000 angle distribution P(OOO) between 
bonded triplets and the structure factor 5'(g), at the low¬ 
est T we have been able to equilibrate and at the optimal 
network density. Previous studies have shown that the 
width of P(OOO) as well as the amplitude of the pre¬ 
peak in S{q) correlate with bond flexibility [s^. Fig. [2] 
shows that the angular distribution is centred around the 
tetrahedral angle but widens on increasing 0max, indicat¬ 
ing the larger number of geometrical arrangements avail¬ 
able for the formation of the network. Simultaneously, 
the larger disorder in the network structure decreases the 
intensity of the S (q) pre-peak, in full agreement with re¬ 
sults based on tetrahedral patchy colloids [s^ . 

To estimate the location of the LL critical point we 
perform grand-canonical Monte Carlo simulations for 
different T to estimate the probability P{N) of find¬ 
ing N particles in the simulated volume V (8 nm^). 
By implementing the successive umbrella sampling tech- 



FIG. 1. Schematic plot of the ST2 water model and of the pro¬ 
posed extension to modulate hydrogen bond flexibility. Solid 
lines indicate the position of the H and LP sites in the rigid 
original ST2 model. The cones have an angular amplitude 
equal to 0max and define the volume limiting the position of 
the same sites in the flexible model (dashed lines). 



FIG. 2. (a) Probability distribution of the OOO angle for 
bonded triplets at p = 0.90 g/cm® for three different values of 
the flexibility. Note the increase of the variance of the distri¬ 
bution on increasing flexibility. Two adjacent water molecules 
are considered bonded if the 00 distance is less than 3.2 A. 
(b) Structure factors S{q) for the same state points, display¬ 
ing the effect of the flexibility on the pre-peak. 
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FIG. 3. Distribution P{N) in the number N of particles pop¬ 
ulating a volume of 8 nm^ at fixed temperature and chemical 
potential /r. This last quantity controls the average density 
and it is selected to provide equal area below the low-density 
and high-density liquid phases. P{N) evolves from a single¬ 
peak to a double peak shape on crossing . The data for 
cos^max = 1 are redrawn from Ref. [l^ . 


nique we have distributed the evaluation of P{N) 
over multiple processors, each of them evaluating the ra¬ 
tio P{N P l)/P{N), for 220 < N < 350. More than 
1000 processors running full time have been dedicated 
for six months to these calculations. Close to a second- 
order critical point, P{N) develops a double peak struc¬ 
ture that becomes more and more pronounced on cooling, 
signalling the coexistence of phases with different den¬ 
sity. During all runs, we have constantly checked that the 
number of crystalline particles (evaluated with the stan¬ 
dard algorithms for detecting ice local structures 23, 3^) 
never exceeded ten nor showed any trend toward growing. 
The results of these calculations for different 0max and T 
are reported in Fig. [31 spanning the T interval over which 
P{N) crosses from a single to a double-peaked function 
with a peak-to-valley ratio around 0.5 (the characteristic 
value assumed by P{N) at the critical temperature [s^) 
down to T where the two peaks are well resolved, sig¬ 
nalling the onset of a clear free energy barrier between 



FIG. 4. Dependence of the liquid-liquid critical point tem¬ 
perature on bond flexibility calculated from free energy 
estimate based on successive umbrella sampling simulations. 
These grand-canonical simulations assume a volume of 8 nm®. 
The dashed red line indicates Tm, the melting T for the liq¬ 
uid to ice Ih/c transformation at the critical pressure. The 
two insets show respectively the critical pressure Pc and the 
critical density pc as a function of 0max. 


the low-density and high-density liquids. The estimated 
location of the critical temperature for the investi¬ 
gated box side {L = 2 nm) as a function of 0max are shown 
in Fig. 01 Consistent with what was previously found for 
the patchy and DNA colloidal models, increasing bond 
flexibility (i.e. increasing flmax) results in lowering . 
Additionally, upon increasing the flexibility, the critical 
pressure decreases and the critical density increases, con¬ 
sistent with the coupling between bond flexibility and lo¬ 
cal density. Indeed, for tetrahedral patchy particles it 
has been shown that the density at zero pressure of the 
fully bonded network decreases with increasing bond di¬ 
rectionality. Similarly, increasing flexibility results in a 
network that is much more easily compressible . Our 
results suggest that the coupling between local density, 
compressibility and flexibility also affects the critical pa¬ 
rameters. 

To properly frame the LL transition in terms of ther¬ 
modynamic stability compared to Ih/c we evaluate the 
free energy of the liquid and of the two ices. To eval¬ 
uate the liquid free energy, we perform thermodynamic 
integration from the ideal gas [28| . To evaluate the crys¬ 
tal free energy we integrate from an Einstein crystal in 
the molecular framework (^ . extending the method to 
account for the flexible arms. For this, we use as a refer¬ 
ence system a thermalised ice lu or Ic configuration with 
at least 20000 particles to average over proton-disorder. 
For each molecule in the reference system we define the 
reference oxygen position (ro), the reference orientation 
of the dipole and HH unit vectors and the reference ori¬ 
entation of the OH and OLP unit vectors, all in the ideal 
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tetrahedral geometry. For each particle, the reference 
Hamiltonian consists of three parts: 


-^trans — 

r-ro|)>^ 

Hiot — 

sin^ pa + (—'] 


\tt J 


4 

.^arms — ^ ^ [1 COS ( 3 ) 

i=l 

where |r — ro| is the distance between the position of the 
oxygen atom in the reference and in the instantaneous 
configuration, a = Inm is a convenient length scale, tpa 
and (j)b are respectively the angle between the reference 
and the instantaneous position of the ideal dipole and HH 
unit vectors and 9i is the angle between the instantaneous 
position of the i patch unit vector and the ideal position 
of that unit vector. In other words, this is the angle we 
limit in our model when specifying 0max- 

Monte Carlo moves, preserving the centre of mass 
position are performed by randomly translating a 
molecule, rotating a single patch (which changes only 
Hurms), or rigidly rotating the water molecule (which 
changes only i?rot)- 

The reference free energy (per particle) of the fully 
constrained system (limit of high A) is (with /3 = l/ksT 
and ks the Boltzmann constant): 

Pfref = /3/trans + /5/rot + 4/3/arm, (4) 


with 


/3/trans — ^ log 


3(Af-l)/2 


V 


/3/rot = - log [v^/4] + - l0g(/3Ar) 


/3/arm = “ log 


(1 - exp[-/3Aa(l - COSdmax)]) 


(1 COS^max)/3AQ 

~ l0g((l - COs6»max)/3Aa) 


( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 


The model free energy / is then calculated following 
the methodology described in Ref. Q. As for the origi¬ 
nal ST2 model, at all temperatures and for all values of 
^max, we find that ice R and Ic have the same free ener¬ 
gies, within our numerical precision (with an uncertainty 
in Pf of ±0.01). From the free energy and the equation of 
state we evaluate the chemical potential = fif+pP/p, 
where P is the pressure and p the number density. The 
main results of the Letter are reported in Fig. [3 showing 
the P dependence of the liquid and ice ^hjc. chemical po¬ 
tential at the LL critical temperature In the case of 
the original ST2 model (dmax = 0°), /3/x of is always 
lower than the liquid one, consistent with the location of 
in the no-mans land. On increasing dmax, the rela¬ 
tive stability of the ^hjc compared to the liquid changes. 
At 0inax — 8°, the liquid chemical potential is slightly 


lower than ice \hlc , while for dmax = 11-5°, at the 
liquid phase has gained a significant stability compared 
to the open crystal lattice. Since the crystal free energy 
is higher than the liquid one (already for dmax — 8° !) 
there is no possibility that the low density liquid phase 
in flexible ST2 will ever convert into the Ih/c structure. 
The low-density liquid phase is, without any ambiguity, 
a phase W itself, definitively disproving the arguments 


Ref. 


ie by itsen 

. [Hill. 


phase diagram, in which we include T, P and 0max, 


In an expanded representation of the 

the 

lines of LL critical points moves with continuity from a 
condition of metastability with respect to Ih and Ic to 
a condition of stability, around 0max = 8°. This con¬ 
tinuity allows us to conclude that the LL critical point 
observed in the original ST2 model must also be gen¬ 
uine. We stress that our study does not aim at providing 
a (better or worse) model for water but to show — with 
an extremely simple modification to the ST2 model — 
that the liquid-liquid transition can be made thermody¬ 
namically stable. For the case of water, the estimated LL 
critical point is definitively located in the region in which 
ice nucleation in bulk is dominant. There, only ingenious 
experiments in the negative pressure region of the phase 
diagram or in the low T glass phases Q can pro¬ 

vide important clues. Still, our proof reinforces the idea 
that the LL transition is a genuine phenomenon in all 
tetravalent systems for which a suitable softness allows 
for network interpenetration and a suitable bond flexi¬ 
bility allows for enhanced stability of the liquid phase(s) 
compared to open crystal structures. 
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